%--------------------------------------------------------------------------
% Return moments under Q
%--------------------------------------------------------------------------

k_vec       = repmat(k',[nS,1]);
k_vec       = k_vec(:);

kb_vec      = repmat(k_bar,[nK,1]);
kb_vec      = kb_vec(:); % [nS*nK_fine,1]

mu_vec      = repmat(mu_XQ,[nK,1]);
mu_vec      = mu_vec(:); % [nS*nK_fine,1]
sig_vec     = repmat(sig_X,[nK,1]);
sig_vec     = sig_vec(:); % [nS*nK_fine,1]
Q_mat       = repmat(Q,[nK,1]); % [nS*nK_fine,nS]

default     = reshape( k' <= k_def ,[nK*nS,1]);
issuance    = reshape( k' >= k_iss ,[nK*nS,1]);

EQ_S        = zeros(nS,nK);
EQ_S2       = zeros(nS,nK);

EQ_Sx        = zeros(nS,nK);
EQ_Sx2       = zeros(nS,nK);

[Cheb_nodes,~] = qnwcheb(nE,-1,1);
a           = -5;
b           = 5;
Eshocks     = (b-a)*(Cheb_nodes'+1)/2 + a;
temp        = pi*(b-a)/(2*nE);
weights     = temp.*sqrt(1-Cheb_nodes'.^2);
probE       = normpdf(Eshocks).*weights;

% kappa transition
g               = exp(mu_vec + sig_vec.*Eshocks);
kp              = k_vec  .* g; %
kp_issuance     = kb_vec .* g; %
kp(issuance,:)  = kp_issuance(issuance,:);


for is=1:nS % future state
    lamS_fit    = griddedInterpolant(k,lamS(is,:)','linear');
    lamS_1      = lamS_fit(kp);
    
    cash_flow   = g.*lamS_1;
    EQ_S        = EQ_S + reshape(Q_mat(:,is).*sum(cash_flow.*probE,2),[nS,nK]); 
    EQ_S2       = EQ_S2 + reshape(Q_mat(:,is).*sum(cash_flow.^2.*probE,2),[nS,nK]); 
    
    lamS_fit    = griddedInterpolant(k,lamSx(is,:)','linear');
    lamSx_1     = lamS_fit(kp);
    
    cash_flow   = g.*lamSx_1;
    EQ_Sx       = EQ_Sx + reshape(Q_mat(:,is).*sum(cash_flow.*probE,2),[nS,nK]); 
    EQ_Sx2      = EQ_Sx2 + reshape(Q_mat(:,is).*sum(cash_flow.^2.*probE,2),[nS,nK]); 
end

EQ_R = EQ_S ./ lamSx;
EQ_R(default) = NaN;

EQ_R2 = EQ_S2 ./ lamSx.^2;
EQ_R2(default) = NaN;

VQ_R = EQ_R2' - EQ_R'.^2;
SQ_R = sqrt( VQ_R );

EQ_Rx = EQ_Sx ./ lamSx;
EQ_Rx(default) = NaN;

EQ_Rx2 = EQ_Sx2 ./ lamSx.^2;
EQ_Rx2(default) = NaN;

VQ_Rx = EQ_Rx2' - EQ_Rx'.^2;
SQ_Rx = sqrt( VQ_Rx );

%--------------------------------------------------------------------------
% Return moments under P
%--------------------------------------------------------------------------

k_vec       = repmat(k',[nS,1]);
k_vec       = k_vec(:);

kb_vec      = repmat(k_bar,[nK,1]);
kb_vec      = kb_vec(:); % [nS*nK_fine,1]

mu_vec      = repmat(mu_X,[nK,1]);
mu_vec      = mu_vec(:); % [nS*nK_fine,1]
sig_vec     = repmat(sig_X,[nK,1]);
sig_vec     = sig_vec(:); % [nS*nK_fine,1]
P_mat       = repmat(P,[nK,1]); % [nS*nK_fine,nS]

default     = reshape( k' <= k_def ,[nK*nS,1]);
issuance    = reshape( k' >= k_iss ,[nK*nS,1]);

EP_S        = zeros(nS,nK);
EP_S2       = zeros(nS,nK);

EP_Sx        = zeros(nS,nK);
EP_Sx2       = zeros(nS,nK);

% kappa transition
g               = exp(mu_vec + sig_vec.*Eshocks);
kp              = k_vec  .* g; %
kp_issuance     = kb_vec .* g; %
kp(issuance,:)  = kp_issuance(issuance,:);

for is=1:nS % future state
    lamS_fit    = griddedInterpolant(k,lamS(is,:)','linear');
    lamS_1      = lamS_fit(kp);
    
    cash_flow   = g.*lamS_1;
    EP_S        = EP_S + reshape(P_mat(:,is).*sum(cash_flow.*probE,2),[nS,nK]); 
    EP_S2       = EP_S2 + reshape(P_mat(:,is).*sum(cash_flow.^2.*probE,2),[nS,nK]); 
    
    lamS_fit    = griddedInterpolant(k,lamSx(is,:)','linear');
    lamSx_1     = lamS_fit(kp);
    
    cash_flow   = g.*lamSx_1;
    EP_Sx       = EP_Sx + reshape(P_mat(:,is).*sum(cash_flow.*probE,2),[nS,nK]); 
    EP_Sx2      = EP_Sx2 + reshape(P_mat(:,is).*sum(cash_flow.^2.*probE,2),[nS,nK]); 
end

EP_R = EP_S ./ lamSx;
EP_R(default) = NaN;

EP_R2 = EP_S2 ./ lamSx.^2;
EP_R2(default) = NaN;

VP_R = EP_R2' - EP_R'.^2;
SP_R = sqrt( VP_R );

EP_Rx = EP_Sx ./ lamSx;
EP_Rx(default) = NaN;

EP_Rx2 = EP_Sx2 ./ lamSx.^2;
EP_Rx2(default) = NaN;

VP_Rx = EP_Rx2' - EP_Rx'.^2;
SP_Rx = sqrt( VP_Rx );

VRP = VQ_R - VP_R;
SRP = SQ_R - SP_R;
